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Pulsed gamma rays have been detected with the Fermi Large Area Telescope (LAT) from more than 20 millisec- 
ond pulsars (MSPs), some of which were discovered in radio observations of bright, unassociated LAT sources. 
We have fit the radio and gamma-ray hght curves of 19 LAT-detected MSPs in the context of geometric, outer- 
magnetospheric emission models assuming the retarded vacuum dipole magnetic field using a Markov chain 
Monte Carlo maximum likelihood technique. We find that, in many cases, the models are able to reproduce the 
observed light curves well and provide constraints on the viewing geometries that are in agreement with those 
from radio polarization measurements. Additionally, for some MSPs we constrain the altitudes of both the 
gamma-ray and radio emission regions. The best-fit magnetic inclination angles are found to cover a broader 
range than those of non-recycled gamma-ray pulsars. 



1. INTRODUCTION 

Millisecond pulsars (MSPs) are thought to be old, 
recycled objects which have reached short spin peri- 
ods 10 ms) via accretion from a binary compan- 
ion Q. MSPs have been established as a population 
of high-energy (HE, > 0.1 GeV) emitters via the de- 
tection of significant pulsed gamma-ray signals from 
more than 20 MSPs at the radio periods using data 
from the Fermi Large Area Telescope (LAT). Addi- 
tionally, more than 30 previously unknown MSPs have 
been found in radio observations of unassociated LAT 
sources, some of which have already been confirmed 
as gamma-ray pulsars [1, [l^ HI] . MSPs can be sorted 
into three sub-classes: those for which the gamma-ray 
peaks lag the main radio component (class I), those for 
which the gamma-ray and radio profile components 
are aligned in phase (class II), and those which lead 
the main radio component (class III). Non- recycled 
gamma-ray pulsars do not show such a diverse distri- 
bution of light curves; thus, studies of MSP gamma- 
ray and radio light curves provide a much broader 
window into the pulsar magnetosphere and emission 
mechanisms. We have generated simulations of dif- 
ferent emission models assuming a vacuum retarded 
dipole magnetic field geometry ^ and attempted to 
reproduce the gamma-ray and radio light curves of 19 
MSPs detected with the Fermi LAT using 2 years of 
data and a maximum likelihood fitting technique. 



2. GAMMA-RAY EMISSION MODELS 

To reproduce the gamma-ray light curves of MSPs 
in class I we have used two-pole caustic (TPC; [6]) 



and outer gap (OG; e.g., [8|) models. Note that we 
take the TPC model to be a geometric realization of 
the slot gap (SG) model [l^ and thus some of the 
parameters of our TPC model differ from those of @ 
as discussed in Section H) 

The TPC/SG model assumes that particles are ac- 
celerated and HE gamma rays are emitted in narrow 
gaps along the surface of last-closed field lines. The 
gaps start at the stellar surface and continue out to 
the light cylinder (defined by the cylindrical radius 
Rlc = c/Q, where is the rotational frequency 
of the pulsar). The OG model assumes that a vac- 
uum accelerating gap forms along the surface of last- 
closed field lines but only above the null-charge surface 
(NCS, defined by the requirement that H. ■ B — 0). 
In the OG model, HE gamma rays are emitted in a 
thin layer interior to the accelerating gap. 

In both models the bulk of the HE emission origi- 
nates at high altitudes above the stellar surface, near 
the light cylinder. The bright, sharp peaks observed 
in HE pulsar light curves [l[ are thought to be emis- 



sion caustics [22| where relativistic aberration, time- 
of-flight delays, and magnetic field line curvature com- 
bine to cancel out phase differences between photons 
emitted at different altitudes. This results in many 
photons arriving at the observer at nearly the same 
phase. 

The TPC and OG models have also been used to 
reproduce the gamma-ray light curves of those MSPs 
in class II; however, given the alignment in phase of 
the gamma-ray and radio peaks for these pulsars we fit 
for the maximum radial extent to which the gamma- 
ray emission is followed, see Sections |4] and [5] for more 
details. We identify these models as altitude-limited 
TPC and OG (alTPC and alOG, respectively) in order 
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to distinguish them from the models used for MSPs 
in class I, for more details on the alTPC/OG models 
see dl. 

The periods (P) of MSPs are found to be very sta- 
ble with small period derivatives (P < 10~^^ s s~^). 
Such relatively low spin-down rates place the major- 
ity of MSPs below the pair-creation death line on a 
P vs. P plot Pulsars below this death line were 
not expected to be capable of screening the acceler- 
ating electric field over most of the open volume and 
thus creating narrow accelerating gaps required in SG 
and OG models. Thus, [ll| created the pair-starved 
polar cap (PSPC) model for HE gamma-ray emission 
in pulsars below the death line. The PSPC model has 
been used to reproduce the gamma-ray light curves of 
MSPs in class HI. In this model the entire open vol- 
ume (from the magnetic dipole axis to the surface of 
last-closed field lines) is available to accelerate parti- 
cles which emit HE gamma rays. 



3. RADIO EMISSION MODELS 

To reproduce the radio light curves of MSPs in 
classes I and HI we have used a single-height, hollow- 
cone beam following the description of [26] and [16]. 
This model assumes that the radio emission origi- 
nates at an altitude which depends on P, the fre- 
quency of emission, and weakly on P [2Q] . For typical 
MSP periods and observing frequencies this altitude 
is < 30% Rlc- This is significantly lower in altitude 
than the bulk of the HE emission and thus leads to the 
non-zero phase offsets between radio and gamma-ray 
light curve features for MSPs in classes I and HI. 

Radio and gamma-ray light curves with features 
which are aligned in phase imply, at least partial, co- 
location of the emission regions. Thus, to reproduce 
the radio light curves of MSPs in class II we used 
alTPC/OG models in which we fit both the minimum 
and maximum radio emission altitudes. This implies 
that the radio emission is also caustic in nature which 
has important implications for the predicted polariza- 
tion properties [28^] . 



4. LIGHT CURVE SIMULATIONS 

We have generated simulations with spin periods of 
1.5, 2.5, 3.5, 4.5, and 5.5 ms for MSPs in classes I and 
HI. For MSPs of class II we have only generated sim- 
ulations with a spin period of 1.5 ms. This is shorter 
than the observed periods of PSRs J2214-^3000 ,25] 
and J1823— 3021A [11| but will, at most, slightly over- 
estimate any predicted off-pulse interval. 

All of our simulations have a resolution of 2° in 
pulse phase and 1° in both magnetic inclination angle 
(a) and observer viewing angle (C). For MSPs in class 



I our simulations have a resolution of 5% of the polar 
cap opening angle for the accelerating and emitting 
gap widths. For MSPs in class II our simulations have 
the same resolution in gap widths and a resolution of 
O.IRlc in emission altitude. 

Our simulation code follows that of 7] with a 
few important modifications. We have included the 
Lorentz transformation of the magnetic field from an 
inertial observer's frame to the frame which co-rotates 
with the star before calculating direction tangent to 
the field line along which a photon is emitted (see 
Appendix B of for more details) as advocated by 
[4I| for self-consistency. Additionally, for the PSPC 
models we have used the same functional form for the 
accelerating field as [13] to calculate the number of 
photons emitted at each step along the field lines as 
opposed to assuming uniform emissivity as is done 
for the other models. For all models the emission is 
never followed beyond a radial distance of 1.2Rlc or a 
cylindrical distance of 0.95Rlc, whichever is reached 
first. In this respect our TPC models differ from those 
originally used by [6] as they only followed emission 
out to a radial distance of 0.95Rlc but not beyond a 
cylindrical distance of 0.75Rlc- 



5. LIKELIHOOD FITTING 

We have developed a Markov chain Monte Carlo 
(MCMC) maximum likelihood technique in order to 
pick the best-fit model parameters. An MCMC proce- 
dure involves taking random steps in parameter space 
and accepting or rejecting new steps based on some 
criteria which only involves the previous step (in our 
case it is the likelihood ratio between the two states) . 
In order to speed up chain convergence and mixing 
we have implemented small- world chain steps [ll| and 
simulated annealing [2lj. We use Poisson likelihood 
for the gamma-ray light curves and a statistic for 
the radio, combine them, and then maximize the joint 
likelihood. 

The formal uncertainty on the radio profiles is much 
smaller, relatively, than that of the gamma-ray light 
curves which drives the likelihood to favor the radio 
fit. In order to balance the relative contributions to 
the likelihood from the radio and gamma-ray data we 
have chosen to use the same uncertainty for each bin 
of the radio light curve which is defined as follows. 
We first calculate the average relative error in the on- 
peak interval of the gamma-ray light curve. Then, we 
multiply that value with the maximum of the radio 
profile. Finally, in the event that we use more bins in 
the radio profile than the gamma-ray light curve we 
multiply the latter uncertainty by the ratio of gamma- 
ray to radio bins. It is important to note that using 
a different radio uncertainty can affect the best-fit ge- 
ometry, in some cases leading to changes in either a 
or C of ~30°. 
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6. RESULTS 

Example fits from each MSP sub-class are shown 
in Fig. [TJ The gamma-ray events for each MSP were 
required to be within 0?8 of the radio position, have 
reconstructed energies from 0.1 to 100 GeV, and have 
zenith angles < 105°. The background levels for the 
gamma-ray light curves were estimated using the LAT 
Science Tool gtsrcprob and spectral results from a 
preliminary version of the 2FGL catalog^ while the 
radio backgrounds are estimated by fitting a constant 
value to the off-pulse intervals. In Fig. [1] the model 
light curves corresponding to TPC and alTPC fits are 
shown in pink, the OG and alOG in green, and the 
PSPC in blue. 




Figure 1: Example light curve fits representing the three 
different MSP sub-classes. Gamma-ray data and model 
light curves are shown in the top panels while the radio 
data and model light curves are shown in the lower 
panels. (Left:) Light curves for PSR J0030+0451 fit with 
OG and TPC models (class I). (Middle:) Light curves for 
PSR J0034-0534 fit with the alOG and alTPC models 
(class II). (Right:) Light curves for PSR J1744-1134 fit 
with the PSPC model (class III). 

When plotting the best-fit geometries (Fig. ^ an 
interesting trend appears. The C values tend to pre- 
fer higher angles near 90°; this is consistent with the 
assumption of a random angular distribution (i.e., 
weighted by sin((^)) of spin axes with respect to the 
Earth line-of-sight. However, the a values seem to 
favor all angles equally in contrast to what has been 
found for non-recycled gamma-ray pulsars [IJ]. This 
may be a manifestation of the recycling process tend- 
ing to align the magnetic and spin axes. 

There are some suggestions in our fits that the 
assumed radio emission altitude of the conal model 
(classes I and III) is too low. One such example is 
PSR J2302-I-4442 for which the gamma-ray model re- 
produces the observed light curve well but the radio 



^ http: / / fermi.gsfc.nasa.gov/ssc / data/ access / 
lat/2yr .catalog/ 
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Figure 2: Best-fit (C,q) values for each MSP (both fits 
are shown for MSPs in classes I and II). TPC fit values 
are shown as filled pink squares, OG as filled green 
triangles, and PSPC as filled blue triangles. The alTPC 
fits are shown as open pink squares and the alOG as 
open green triangles. 



model can not produce widely enough separated peaks 
at the best-fit geometry Q. Fig. [3] demonstrates how, 
for the same geometry which gives a good gamma- 
ray fit, increasing the radio emission altitude does, in 
fact, lead to a better fit to the radio data. Note that 
the radio altitude was increased in the simulations by 
decreasing the emission frequency. 
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Figure 3: Radio data and models with different emission 
altitudes for PSR J2302+4442. 



7. FUTURE 

The simulations described in Section 2] have rather 
coarse resolution in the gap width and altitude param- 
eters and we have found that a finer gridding is needed 
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before we can look for meaningful trends in these val- 
ues. In some of the MSPs we have analyzed the single, 
hollow-cone beam model is not correct as there is ev- 
idence for either a core component or multiple cones 
from radio polarization data. Additionally, as noted 
in Section |6l we find indications that the radio emis- 
sion should originate at higher altitudes in the MSP 
magnetospheres. Therefore, in addition to increasing 
the resolution of our simulations we plan to produce 
more complex radio models and explore changing the 
emission altitude in order to more closely match the 
radio profiles. We will also compare the predicted po- 
larization angle swings with polarimetric data to guide 
further model refinement and serve as an additional, 
observational test. 

The magnetosphere of a pulsar should be filled, to 
some extent, with charges (e.g., jl^); thus, the vac- 
uum solution for the magnetic field can not exactly 
match reality. With that in mind, we plan to apply the 
same fitting technique to simulations using magnetic 
field geometries from magneto-hydrodynamic simula- 
tions of a pulsar magnetosphere under force-free as- 
sumptions (e.g., and with finite conductivity [Tsj . 
By comparing our predicted gamma-ray and radio 
light curves in different magnetic field geometries with 
those from observations, we can constrain the struc- 
ture of a real pulsar magnetosphere. 
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